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1 . I nt  roduct ion 


Finite  difference  methods  specify  the  dependent  variables  at 
certain  grid  points  in  space  and  time, and  the  derivatives  in  the  equations 
are  evaluated  using  Taylor  series  approximations.  The  Galerkin  procedure, 
which  will  be  treated  in  this  chapter,  represents  the  dependent  variables  with 
a sum  of  functions  which  have  a prescribed  spatial  structure.  The  coefficient 
associated  with  each  function  is  normally  a function  of  time.  This  procedure 
transforms  a partial  differential  equation  into  a set  of  ordinary  differential 
equations  for  the  coefficients.  These  equations  are  usually  solved  with 
finite  differences  in  time.  The  two  most  useful  Galerkin  methods  are  the 
spectral  method  and  the  finite  element  method.  The  spectral  method,  which 
employs  orthogonal  functions,  has  been  used  in  meteorological  problems  for 
a number  of  years.  The  finite  element  method  employs  functions  which  are 
zero  except  in  a limited  region  where  they  are  low  order  polynomials.  This 
method,  which  was.  developed  in  engineering,  has  only  recently  been  introduced 
into  meteorology  and  oceanography. 

The  Galerkin  procedure  can  be  illustrated  with  the  following  equation: 

(u)  - f (x)  (1) 


where  £ is  a differential  operator,  u is  the  dependent  variable  and  f(x) 
is  a specified  forcing  function.  Suppose  that  (l)  is  to  be  solved  in  the 
domain  a x _<  b and  that  appropriate  boundary  conditions  are  provided.  Con- 
sider a series  of  linearly  independent  functions  (^(x)  which  will  be  called 
basis  functions.  The  next  step  is  to  expand  u(x)  into  a series  as  follows: 


N 

UjfPj(x)  , 

J-l 


(2) 
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where  is  the  coefficient  for  jth  basis  function.  The  error  in  satis- 

fying the  differential  equation  (1)  with  the  N terms  of  the  sum  (2)  is 


*n  uj<pj)  ■ f(x)  ■ 


The  Galerkin  procedure  requires  that  the  error  be  orthogonal  to  each  basis 
function  in  the  following  sense: 


b 

/ 


qpjdx  = 0 , 1*1, . . . ,N  . 


The  final  form  is  obtained  by  substituting  (3)  into  (4) 


b N b 

J *1*'^  ' J <pi 


f(x)dx  = 0 , i=l N 


This  reduces  the  problem  to  N algebraic  equations  which  relate  the  unknown 
coefficients  to  the  "transforms"  of  the  forcing  function.  This  procedure 

is  quite  general  and  can  be  applied  to  more  dependent  and  independent 
variables. 

2 Example  with  Spectral  and  Finite  Element  Methods 

Now  the  spectral  method  and  the  finite  element  method  will  be  applied 
to  the  following  simple  form  of  (1): 


f(x)  , o < x < TT  . 


The  boundary  conditions  are 


u(o)  ■ u(tt)  = 0 . 


For  the  spectral  method  the  following  basis  functions  are  appropriate: 
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sin  Jx  , J-l N . 


(8) 


These  functions  are  orthogonal  on  the  Interval  0 < x ti  and  they  satisfy 
the  boundary  conditions  (7).  With  these  basis  functions 

N V 

£.(Y  ‘“j'j”  ■ Yu  (~j2,v'.i  • 


j-i 


and  ( 5)  becomes 


__  \ ^ i *•..  /*  j Cp4 d x * / (p.  f ( x ) d x , i*l,  • • • |N  • 

k SJ  J 


(9) 


The  product  of  the  basis  functions  can  be  written 


ff 

/ 


slnlx  slnjx  dx 


v 

*/ 


[ cos  ( l- j ) x-cos  ( 1+J  ) x ] dx  - ( tt/ 2)  6 ^ , 


(10) 


where  6^  Is  the  Kronecker  delta  which  satisfies  6 


1 if  i - J and 


- 0 if  i t J . Equation  ( 10)  is  merely  the  orthogonality  condition 
which  arises  since  the  Integral  vanishes  except  when  i - j . With  the  use 
of  (10),  the  solution  to  (9)  becomes 

n 

2 


n 


/ 


f dx  . 


(11) 


Each  coefficient  is  proportional  to  the  finite  Fourier  transform  of  the  forc- 
ing term.  In  this  example  both  the  error  in  the  solution  and  the  error  in 
the  differential  equation  are  orthogonal  to  the  basis  functions.  This  is 
because  of :«i)  is  proportional  to  (j\  so  that  if  the  error  is  orthogonal  to 
JC((?  ) it  will  also  be  orthogonal  to  CT  . This  will  also  be  true  when  certain 
other  linear  equations  are  treated  with  the  spectral  method,  but  it  will  not 


generally  be  true  with  nonlinear  equations 
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Now  consider  the  same  differential  equation  (6)  with  the  finite  ele- 
ment method.  Divide  the  interval  0 < x < it  into  N+l  segments  such  that 
(N+l)  Ax  - tt  . The  basis  functions  are  chosen  to  be  tent  shaped  piecewise 
linear  functions  which  are  also  called  chapeau  functions,  as  shown  in 
Figure  1.  As  can  be  seen  from  the 


x -*■ 

Fig.  1.  Piecewise  linear  basis  function. 


Figure,  tPj (x)  has  a maximum  of  1 at  x = jAx  , which  is  called  the  nodal 
point.  The  basis  function  decreases  linearly  to  zero  at  x = (j-l)Ax  , and 
it  is  zero  everywhere  else.  Mathematically  cp^  (x)  is  defined  as  follows: 

0 , x > (j+l)Ax  or  x < (j-l)Ax  I 

cpj  (x)  = (x- (j-l)Ax) / Ax  , (j-l)Ax  _<  x £ jAxl 

( (jtl)Ax-x) /Ax  , jAx  £ x £ (j+3)Ax| 


Note  that  the  coefficient  u^  is  actually  the  value  of  the  function  at 


x * jAx  since  cpj  (jAx)  * 1 and  cp^(jAx)  = 0 for  i f*  j . These  elements 
are  not  quite  orthogonal,  but  only  adjacent  elements  interact.  The  boundary 
conditions  (7)  are  automatically  satisfied  although  this  is  not  necessary 
in  many  cases  with  finite  elements. 


Equation  (5)  now  becomes 


This  form  of  the  equation  is  not  appropriate  because  it  involves  a second 

derivative  of  the  basis  function  which  is  only  piecewise  linear.  However, 

this  problem  can  be  avoided  by  integrating  the  first  term  by  parts  as  follows: 
IT  IT 


x)dx  = 0 . 


J = 1 o 


The  first  term  vanishes  because  all  of  the  op's  are  zero  at  x = 0,  TT  . 


The  Galerkin  equation  now  becomes 
IT  TT 


-i'J 


d<^  dcp 


dx  dx 


dx 


f (x)dx  , i-1, . . . ,N  . 


(13) 


j-1  o 

Note  that  differentiating  (12)  gives: 

0 , x >(j+l)Ax 


dcp 

dx 


or  x < (j-l)Ax 


( 1/Ax  , (j-l)Ax  £ x £ JAx 
-1/Ax  , jAx  £ x £ (j+l)Ax 


(14) 


The  left  hand  side  of  (13)  is  easily  evaluated  since  only  3 terms  in  the 


sum  are  different  from  zero: 
TT 
N 


u,  ,Ax  - 2u.Ax  + u^.-Ax 


-E'J 


dx 


(15) 


Ax 


j=l  o 

We  right  hand  integral  in  (13)  may  be  evaluated  if  f (x)  is  approximated 
in  terms  of  the  basis  functions: 

N 


f(x>  = fjcpj  * 
J=1 


(16) 


so  that  the  integral  becomes 
TT  IT 


TT  TT  N (i+l)Ax 

/N  N /% 

cp1f(x)dx  = E £i  / E‘i 

i ^ (i--i)A* 
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If  £ =»  x - l Ax  Is  introduced  the  integral  can  be  expanded  into  three  in- 


tegrals: 


Ax 


/rf  x.  r f ^(f,+Ax)  ,r  r (C+Ax)2jc.ic  / f,(Ax-f,) 

tp*f<x)dx  ■ -£*-i  J r/  '1  J ’ J 


dr. 


-Ax 


-Ax 


(17) 

When  these  terms  have  been  evaluated,  ( 17)  and  (15)  can  be  substituted 
into  (13)  which  gives 


“i+l  “ 2ui  + Ui-1 


fi+l  + 4fi  + fi-l 


Ax 


(18) 


This  equation  applies  for  2 _<  i <_  N-l  and  the  equations  for  i = 1 and 
i = N are  obtained  by  removing  any  terms  in  1 = 0 or  i = N+l  . Equation 
(18)  may  be  solved  by  Gaussian  elimination 

Since  each  coefficient  in  this  finite  element  expansion  represents  the 
solution  at  a certain  point  in  space,  it  is  convenient  to  compare  (18) 
with  finite  difference  forms  of  (6).  The  centered  difference  form  (A) 
of  this  equation  is 


UiH  ~ 2ui  * Ui-1 
Ax2 


I ’ 


(19) 


where  ui  = u(iAx)  . The  finite  element  equation  (18)  and  the  finite 
difference  equation  (19)  are  the  same,  except  that  the  forcing  term  in 

(18)  appears  in  a weighted  average.  When  these  equations  are  solved  with 
a f(x),  which  is  sinusoidal,  the  finite  element  form  is  considerably  more 
accurate  for  the  shorter  wavelengths. 

In  this  example  it  appears  that  the  spectral  method  is  superior  because 


the  solution  error  is  actually  orthogonal  to  the  basis  functions.  This  is 


not  generally  true  with  the  finite  element  method  because 
L(u  ) depends  on  u^^  , and  u1+1  . Each  Increase  In 

N will  normally  change  all  of  the  solutions  , whereas  with  the  spectral 

method  the  original  N amplitudes  are  not  changed  because  they  are  already 
exact.  However  If  the  variation  of  f should  require  fine  resolution  in 
only  a small  area,  the  finite  element  method  can  easily  be  applied  by  letting 
Ax  vary.  In  this  case  the  spectral  method  would  require  more  elements  be- 
cause its  spatial  resolution  is  uniform.  It  is  also  clear  that  the  finite 
element  method  can  be  used  to  design  better  finite  difference  equations. 

3 Time  Dependence 

In  the  previous  sections  the  Galerkin  procedure  has  been  applied  to  one- 
dimensional equations  which  are  independent  of  time.  The  treatment  of  time 
variation  is  important  for  most  meteorological  prediction  problems.  Consider 
the  following  simplified  equation: 


% (u) 


= 0 , 


(20) 


where  the  operator  ^ may  be  nonlinear.  Expand  u(x,t)  into  a series  as 
follows: 


u(x,t) 


N 

^ 1 u (t)cp'(x)  , 

J-l  J J 


(21) 


where  the  coefficients  Uj(t)  are  functions  of  time  and  the  basis  functions 
(x)  are  functions  of  x . Usually  the  Galerkin  procedure  is  not  applied 
to  the  time  dependence  because  it  is  more  convenient  to  ue  finite  differences 
in  time. 

The  Galerkin  form  of  (20)  is  obtained  by  substituting  (21)  into 
(20),  multiplying  by  cp^Cx)  and  integrating  over  the  domain  as  follows: 


±*J  f 

^>u  V <W,X  v 


. ‘»i 

I - l a 


<p4  . *>\  >‘i*  - o , t-i, 


Tl» I *•  privf.iM  jilvea  N roup W>il  ordinary  differential  equations  lu  the  co- 
«*  t' t It'  lenta  . This  art  ran  he  solved  l>v  I nt  roduc  I n>t  finite  differ- 

ences tn  time. 

The  Importance  of  enerjtv  consor  v lnfc  l lotto  difference  schemes  I'  wall 
known.  Tho  Galerkln  method  I rails  naturally  to  enorqy  conservation  In  equa- 
tion* with  quadratic  eneiqy  variants.  To  show  this,  multiply  ( 2 0 ) hy  ti 
and  Integrate  with  respect  to  x : 


h 

J * J 


u X(u)dx  . 


(«) 


Kor  an  energy  conserving  system,  tho  oporatoi  must  satisfy  tho  condition 
b 

/ u £(u)dx  - II  . (24) 

*a 

whoro  u la  any  function  which  satisfies  the  boundary  conditions.  In  this 
case  (2  1)  becomes  , 


d 

dt 


/' 


d X - l)  , 


(2S) 


which  shows  the  energy  conservat Ion  tor  the  exact  equation.  To  demonstrate 
that  the  same  result  holds  for  the  finite  sum  (21),  multiply  the  ith 
equat Ion  of  12?)  bv  u ^ and  sum  from  l-l  to  l-N  : 


b N N b N N 

j <E^vL(^uiVax " '/  u 


V.)  • (2b) 


IS 


The  integral  on  the  right  vanishes  from  ( 24>  since  the  function  given  by 

21)  satifies  the  boundary  conditions.  Therefore  (26)  can  be  written 
b N 

ulCPl)2/2  dx  - 0 , (27) 

a 

which  expresses  the  energy  conservation  for  the  Galerkin  approximation  to  the 
spatial  variation.  As  with  finite  difference  equations  the  actual  degree  of 
energy  conservation  will  depend  on  the  time  differencing  which  is  used  in 
(22). 

4 Barotropic  Vorticity  Equation  with  Fourier  Basis  Functions 

In  this  section  the  spectral  method  will  be  applied  to  the  barotropic 
vorticity  equation  on  the  beta  plane.  Fourier  basis  functions  are  appropri- 
ate for  the  beta  plane  when  the  fields  are  periodic  in  x and  y . The 
development  of  this  section  closely  follows  Lorenz  (1960).  The  barotropic 
vorticity  equation  may  be  written: 

V2^  + k x • V(V2ii»)  + 8 9ii»/3x  = 0 , (28) 

where  \p  is  the  streamfunction.  Suppose  that  the  fields  are  periodic  in 
both  x and  y so  that 

iHx  + 2n/k;y  + 2TT/S.,t)  - <Kx,y,t)  • (29) 

With  the  beta  plane  geometry  and  the  periodicity  condition,  the  appropriate 
orthogonal  basis  functions  are  of  the  form: 

<p  <x,y)  - .*<•>■»♦"*»>  . (30) 

mn 

These  functions  are  eigensolutions  of  the  equation: 

V2cp  + bcp  - 0 (31) 


» 

i 

I 
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where  the  eigenvalues  are  given  by 


2 2 2 2 

b - (m  k + n t ) . (32) 

The  streamfunct ion  can  be  expanded  In  terms  of  this  basis  functions  as 
follows : 

£ cnn(0  e‘<»k^«y)  _ 

la  n 


In  order  for  ip  to  be  real  the  coefficients  must  satisfy  the  condition 


* 

C - C 
mn  -m-n 


where  ( ) Indicates  the  complex  conjugation.  This  can  be  shown  by  con- 

sidering only  the  m,n  and  -m,-n  . It  Is  convenient  to  Introduce  the  wave 
number  vector  M - mkt  + ntj^  and  the  radius  vector  R - xt  + yj  . The 
expansion  for  can  now  be  written 


With  the  use  of 


<Mx,y,t) 

(31)  and 


1M*  R 
e 


the  vorticlty  can  be  written 


(33) 


T.  (M-M)  Cnj(t)  e1***  . 
M 


(34) 


The  quantities  which  are  required  in  the  nonlinear  term  in  (?8)  may  be 
written: 

V*  - 1H  C>  e1H** 

H 

V(V2i^)  - IL(L-L)  C+  elL‘®  . (35) 

L 
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The  wavenumber  vectors  H and  L are  Introduced  because  the  suras  must  be 
multiplied  together  and  rearranged. 

Now  substitute  the  various  sums  [(  33),  (3M  and  (35)1  into  (28) 

which  gives: 


(L‘L)^r  C ? e1^'*  Xj  (L*L)k*HxL  C.tC>  e1(“+L)'R 


'dt  L 


H L 


L H 


+ lmg 


£ 


1L*$ 

C*  e - 0 

la 


(36) 


The  Galerkln  method  fos  this  equation  Is  similar  to  the  method  used  In 

(22),  except  that  the  equation  must  be  multiplied  by  the  complex  conjugate 

of  the  basis  function  since  the  basis  function  Is  complex.  To  carry  out  this 

1M*$ 


process  multiply  (36)  by  e 
as  follows: 

2ir/k  2 it/I 


and  Integrate  over  the  periodic  domain 


II 


"E  <£•£>  57  H «1<l'M,‘R  + i«»E  C* 


l (L-M) • R 


/ , / , (L*L)k*HxL  C>Cj-  e 


-*  ->  -4  -f 

i (H+L-M)  • R 


->  > 

L H 


dydx  ■ 0 , 


for  each  M in  the  original  sum  (33).  Each  Integral  of  the  exponential 
function  will  vanish  except  when  the  exponent  is  zero.  This  leads  to  the 
following  equation  for  each  M t 


dCj*  lmkfiCj* 
dt“  + 


sa*  E 

M*M 


(M-H) • (M-H)k*HxM 
M*M 


* C*  - 0 . 


SS-H  CH 


(37) 
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In  the  first  two  terms  the  contribution  occurs  for  L - M and  in  the  last 
term  for  L ■ M - ft  . 

Equation  (3)  represents  N ordinary  differential  equations,  where 
N is  the  number  of  terms  in  the  sum  (33).  The  last  term  in  the  equation 
gives  the  interaction  between  different  waves  which  comes  from  the  nonlinear 
advectlon  term  in  (28).  In  particular  wave  M is  affected  by  the  inter- 

action of  waves  H and  M - H . When  the  last  term  Is  dropped,  (3) 
becomes  a set  of  linear,  uncoupled  equations  which  can  be  solved  to  give 
the  Rossby  wave  solution. 

In  section  (3)  it  was  pointed  out  that  the  Galerkin  procedure  pre- 
serves energy  type  univariants  which  arise  from  quadratic  nonlinearities  in 
the  original  equations.  Equation  (28)  conserves  both  kinetic  energy  and 
mean  square  vorticity  or  enstrophy.  The  kinetic  energy  for  the  region  can 
be  written: 

2:r/k  2ir/£  2n/k  2v/l 

K • f f w*  - ill  1^ "•BchcS 

J I H M 

o o J 

where  the  V 4>  product  was  obtained  from  (35)  with  different  summations. 

The  integral  on  the  right  is  nonzero  only  when  H - -M  so  that  the  energy 
can  be  written 

M-MCfiC-ft  “ \ Z^MlSil2  • (38) 

M M 

* 

where  the  condition  C J ■ C J has  been  used  in  the  last  step. 

The  energy  form  in  (38)  is  conserved  (dK/dt  ■ 0)  by  both  the  original 
vorticity  equation  (28)  and  the  spectral  form  (37).  The  conservation  for 
(28)  is  easily  demonstrated,  and  the  conservation  for  (38)  follows 
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from  the  development  in  section  (3).  An  equation  for  the  rate  of  change 
of  energy  In  wave  M can  be  obtained  by  differentiating  C^C  * with  respect 
to  t and  by  using  (37).  The  resulting  equation  shows  that  the  energy 
in  wave  M changes  in  proportion  C>  times  the  amplitudes  of  pairs  of 
interacting  waves.  Thus  if  is  maintained  at  zero,  the  energy  flow  out  of 

the  other  waves  to  it  must  be  zero.  This  shows  in  another  way  that  energy 
will  be  conserved  in  any  set  of  waves  that  might  be  selected  for  sum  ( 33). 
Since  interactions  outside  of  this  set  are  neglected,  aliasing  cannot  occur 
in  a spectral  model.  This  automatically  eliminates  the  nonlinear  computa- 
tional instability  which  occurs  with  finite  difference  equations. 

The  set  of  ordinary  differential  equations  (37)  can  be  integrated 
numerically  with  one  of  the  standard  schemes.  In  fact  Baer  and 
Platzman  (1961)  noted  that  the  linear  terms  in  (37)  can  be  treated  exactly 
so  that  the  only  time  differencing  errors  comes  from  the  nonlinear  terms. 

It  is  clear  that  the  spectral  method  is  much  more  accurate  than  most 
finite  difference  methods  for  the  same  number  of  degrees  of  freedom.  In  par- 
ticular, linear  advection  that  was  examined  is  treated  exactly 
by  the  spectral  method  provided  that  the  initial  field  in  resolved.  Finite 
difference  methods  experience  false  dispersion  since  the  short  waves  move 
too  slowly.  The  spectral  method  has  no  aliasing  because  interactions  in- 
volving shorter  waves  outside  of  the  truncated  set  are  excluded.  On  the  other 
hand,  the  finite  differencing  falsely  reflects  interactions  with  shorter  waves 
back  onto  longer  waves.  With  the  Arakawa  Jacobian  finite  difference  forms, 
this  aliasing  does  not  produce  spurious  energy,  but  it 

does  cause  phase  errors  in  the  interacting  waves.  In  spectral  models  the 
most  important  error  involves  the  neglect  of  interactions  with  wave  components 
which  are  outside  of  the  original  set.  The  neglect  of  these  interactions 
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causes  an  error  in  the  waves  which  are  represented  by  the  basis  functions. 

Thus  although  the  error  in  the  original  equation  is  orthogonal  to  the  basis 
functions,  the  error  in  the  solution  will  occur  in  the  scales  described  by 
the  basis  functions. 

When  the  spectral  method  is  applied  to  a vorticity  equation  such  as 
(28),  a Poisson  equation  for  3iJ</3t  does  not  have  to  be  solved  since  the 
basis  functions  are  eigensolut ions  of  (31).  The  Poisson  equation  must  be 
solved  at  each  time  step  with  finite  difference  methods.  The  biggest  draw- 
back to  this  form  of  the  spectral  equations  is  in  calculating  the  nonlinear 
tern  which  appears  as  the  sura  in  (37).  The  coefficient  preceding 

‘C:J  it*  called  the  interaction  coefficient  and  it  is  usually  computed 
Just  once  and  stored  for  use  during  the  integration  of  the  equation.  The 

problem  is  that  if  there  are  N degrees  of  freedom  the  number  of  operations 

2 

needed  to  compute  the  nonlinear  terms  goes  as  N~  for  this  spectral  model 
as  compared  with  N for  most  finite  difference  methods.  Thus  for  high 
resolution  (large  N),  this  form  of  the  spectral  method  requires  relatively 
larger  computer  time  than  finite  difference  methods.  In  a later  section  a 
method  which  avoids  this  problem  will  be  presented.  However  the  present 
method  is  very  convenient  for  low-order  models.  Loren*  (1*160)  obtained  some 
very  interesting  nonlinear  solutions  with  a 3-component  system.  It  can  be 
seen  from  (37)  that  at  least  3 waves  are  required  for  nonlinear  interaction. 
S Barotroplc  Vorticity  Kquat ion  with  Spherical  Harmonics 

In  this  section  the  spectral  equations  will  be  formulated  for  barotropic 
mot  ion  on  the  sphere.  The  barotropic  vorticity  equation  in  spherical  co- 


ordinates can  be  written: 


3u  3T"  ~ IX 


3VJiK  211  3i(> 
3y  + 2 3X 


0 . 


(39) 


l 


where 


V2 


+ 


1 

i-u* 


(40) 


In  these  equations  X Is  the  longitude  and  u * slnq>  , where  q>  Is  the 
latitude.  The  spectral  method  was  first  applied  in  spherical  coordinates 
by  Sllberman  (1954)  and  the  development  of  this  section  follows  Platzraan 
(I960). 


The  appropriate  orthogonal  basis  functions  are 


Y 

m,n 


<U.  X) 


P 

ra,n 


(y)  e 


lmX 


(41) 


where  P denotes  associated  Legendre  functions  of  the  first  kind  which 

m,n 

are  defined  by 


P 

m,n 


(VO 


,(2nfl)(n-m)l  1/2  ( l-y2)m/ 2 dnHa  2 ..n 
(n+m)!  1 2nnJ  dyn*n  ^ 


(42) 


These  basis  functions  are  spherical  harmonics  which  satisfy  the  equation 


VJ  Y + b Y - 0 , 
m,n  m,n 


(43) 


where,  the  eigenvalues  are  given  by 


b - n(n+l)/a 


(44) 


Here  |ra|  is  the  planetary  wave  number  and  n-|m|  is  the  number  of  zeros 
between  the  poles.  Also  n must  be  greater  than  or  equal  to  |m|  . These 
basis  functions  are  orthogonal  so  that 

2ir  .1  /I  for  (m',n')  - (m,n) 


-L  f f Y Y* 

4tt  / / m,n  ra',n' 


d v 1 d X 


(45) 


0 for  (m',n')  i (m,n) 
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The  streamfunction  can  be  expanded  as  follows: 


M | m | + J 

>P(x,y,t)  - a 2 ^ ^ 

ra=-M  n*=|m| 


ip 

m,n 


(t) 


Y 

m,n 


(*>U)  . 


(46) 


Since  \p  must  be  real  ip  must  satisfy 

m,n  J 


ip 

-m,n 


=s 


(-1) 


* 

>P 

m,n 


(47) 


This  condition  was  'derived  with  the  use  of  the  relation  P * (-1)™  P 

-m,n  m,n 

The  coefficients  <p^  ^ can  be  obtained  from  the  inverse  transform: 


ip 

n,m 


(t) 


<P( A.u, t) 


Y*  dpdX  . 
mtn 


(48) 


The  vorticity  has  the  following  expansion: 


V2i p 


ip 

m,n 


(t) 


Y 

m,n 


(49) 


which  follows  from  (43)  and  (44). 

The  Galerkin  method  is  applied  by  substituting  (48)  and  (49)  into 
(39),  multiplying  by  Y^  ^ and  integrating  with  respect  to  y and  X . 
When  the  conditions  (45)  are  employed  this  equation  reduces  to: 


dip 


m,_n 


2ftm  1 


dt 


The  nonlinear  terms  F 


* 


m,n 
M ! m j + J 


n(n+l)  *m,n  n(n+l)  m,n 
may  be  written 


(50) 


M 


|m2|+J 


m,n 


X!j  , X^  l^m  ,n  ^m  ,n  Mm.n^.n^iiyn, 

m~--M  n^jraj  m2=-M  n2=|m0  ] 1 *■ 


> , 


(51) 
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where  the  Interaction  coefficients  are  given  by 


1 

( n ^ (tij+1  )-n  , (n ,+  l)  ) I 
dP  -1 

ral  ,ni 

n 1 1 ^ . f 


dP 

P (m,P 
m,n  l m^,n^  d\i 


^ ) 


L(m,n;m  ,n ,;m, ,n.)  * ( - m,P  — — i — ±ldu  for  m - m.+  .t« _ 

l l 2 2 I 2 m,.,n~  dp  1 2 


for  m 4 ra^+  m., 


(52) 


In  obtaining  this  result  the  subscripts  1 and  2 were  used  for  expansions 

(48)  and  (49),  respectively.  This  form  for  the  interaction  coefficients 

comes  from  the  fact  that  F changes  sign  when  U>  and  are 

m,n  ral,ni  m7,n2 

Interchanged. 

Equation  (50)  has  the  same  form  as  the  prediction  equation  (37)  for 
the  Fourier  basis  function.  However  the  spherical  coordinate  equation  has 
more  complicated  Interaction  coefficients  because  of  the  integral  involving 
the  Legendre  functions.  It  can  be  shown  by  the  same  method  as  before  that 
energy  is  conserved,  and  Platzman  (I960)  has  also  shown  that  mean  square 
vorticity  is  conserved.  The  spectral  method  applied  to  spherical  (global) 
prediction  has  the  advantage  that  there  are  no  singular  points  whereas  singular 
points  often  cause  problems  with  finite  difference  models.  The  only  major 
disadvantage  in  solving  ( 50)  is  in  the  large  number  of  terms  which  come  from 
the  nonlinear  terms.  This  problem  will  be  treated  in  the  next  section. 


6.  Transform  Method 

In  this  section  a new  method  for  handling  the  nonlinear  terms  In  (50) 
will  be  presented  which  avoids  the  use  of  interaction  coefficients  (see  (51) 
and  (52)).  This  method  was  formulated  independently  by  Orszag  (1970)  and 
Ellasen,  Machenhauer  and  Rasmussen  (1970),  and  it  has  been  reviewed  by  Bourke, 
McAvaney , Puri  and  Thurling  (1977).  The  problem  with  the  interaction  coeffi- 
cient method  for  computing  nonlinear  terms  is  that  it  requires  multiplication 
of  two  series  (together)  which  is  very  time  consuming.  The  transform  method 
sums  the  series  at  certain  spatial  grid  points  and  these  fields  are  multiplied 


together  at  each  point  to  form  the  nonlinear  terms.  Then  the  nonlinear  terms 
must  be  transformed  back  to  spectral  space.  The  usefulness  of  this  process 
is  enhanced  by  the  existence  of  efficient  transform  methods.  In  spherical  co- 
ordinates the  fast  Fourier  transform  is  used  in  longitude  and  the  Legendre  inte- 
grals in  latitude  are  evaluated  by  Gaussian  quadrature.  This  method  is  far 
superior  to  the  interaction  coefficient  method  for  the  sphere. 

The  nonlinear  terms  which  must  be  transformed  mav  be  rewritten  as  follows: 


rfl.  - 1 3V24»  3^  3V2iK  _ I , 3 9 ,3<J>  n2.., 

F(U,A)  2[9 U^T  9A~HrJ  2l9A(9vi  v 9u(3A  V 


It  is  now  convenient  to  define  the  following  quantities  which  are  the  A and 
Cp  velocity  components  multiplied  by  coscp: 

U , v = ! |$  . (54) 

a 9y  a 9A 


When  these  velocities  are  introduced  into  (53)  it  can  be  written  as  follows: 


F(y,\)  - - ^[— ^-^(uv2^)  + ^(W2U>)  1 • 


(55) 


The  velocity  components  (54)  can  be  computed  from  (46)  at  longitude- 
latitude  grid  points,  and  the  vorticity  can  be  obtained  at  the  same  points 
using  49).  The  details  of  the  process  will  be  given  later.  The  products 
UVJii»  and  VVJi{/  can  be  calculated  at  each  grid  point  and  the  resulting  prod- 
ucts can  be  Fourier  analyzed  in  X to  give  the  following  relations: 

m-M 

UV‘>  ■ a V A (y)  elmX  , 

/ a ra 


ra~-M 


(56) 


VV24> 
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E 

m=-M 


B ly)  e 
ra 


ImA 


The  transform  of  F(y,X)  is  given  by 


F 

m,n 


-imA 

e 


P F(y,X)  dydX  . 
ra,  n 


(57) 


The  X integration  in  (57)  can  be  carried  out  by  substituting  (56)  into 
(55)  and  by  inserting  the  result  in  (57),  which  gives: 


F 

m.n 


1 

2 


dB 

A P + -r-S  P ] du  . 
ra  m.n  dy  ra,n 


The  second  terra  can  be  integrated  by  parts  which  gives 
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m dy 
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(58) 
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where  the  condition  B = 0 at  y = ± 1 was  used  to  simplify  the  integral. 

m 

This  condition  follows  since  V is  equal  to  the  actual  velocity  times  coscp. 

The  form  of  F given  by  (58)  is  superior  to  the  earlier  form  because 
m,n 

only  the  known  function  P is  differentiated. 

ra,n 

The  integrand  in  (58)  is  a polynomial  in  y and  the  integral  can  be 

evaluated  following  Eliasen  et  al.  (1970)  by  the  Gaussian  quadrature  formula. 

If  the  integrand  is  denoted  by  Q(y)  , the  formula  gives  the  following 

expression  for  F : 

m,n 


Q0ik)  . 


(59) 


k=l 


In  (59)  the  summation  is  carried  over  K values  of  y^  , where  the  y^'s 

(K) 

are  roots  of  the  Legendre  polynomial  ^ and  are  the  corresponding 

Gauss  coefficients.  The  formula  is  exact  for  any  polynomial  of  degree  smaller 
than  or  equal  to  2K-1  (see  ).  Thus  apart  from 

roundoff  errors,  no  approximation  is  introduced  by  computing  the  integral  when 
a sufficiently  high  value  of  K is  used.  The  maximum  degree  of  Q(y)  can  be 
most  easily  obtained  from  (32). 

Before  discussing  this  process  for  treating  the  nonlinear  terms  in  more 

detail,"  it  is  necessary  to  determine  the  relation  between  J and  M which  must 

be  defined  in  the  sum  (46).  In  rhcr~boidal  truncation  J = M , so  that  each 

latitudinal  mode  has  the  same  numbex  of  waves  in  longitude.  With  triangular 

truncation  J = M - |m|  so  all  basis  functions  which  have  the  same  scale 
2 

n(n+l)/a  , are  either  retained  or  dropped.  Thus  the  mole  with  the  smallest 
latitudinal  scale  has  the  largest  longitudinal  scale.  Most  meteorological 
models  use  the  rhomboidal  truncation  in  part  because  it  gives  better  longitudinal 
resolution.  In  the  remainder  of  this  development,  the  rhomboidal  truncation 
will  be  used. 
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l»  order  to  construct  the  fields  (56)  It  is  necessary  to  obtain  l) 
and  V from  ij>  . First  expand  l)  and  V Into  these  sums: 

M I ml  +M+- 1 


ni“-M  n=*  j m j 

M | m | fM 

v • V V 

in** -M  n-  | m| 


l)  Y 
m ,n  in , n 


V Y 
m , n ra , n 


(60) 


The  following  relations  will  be  useful  In  evaluating  (54): 
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where  D ’ I (n  -m  )/(4n  -1)1  . The  final  expressions  for  U and 

m , n in , n 

V can  be  obtained  bv  substituting  (46)  and  (60)  Into  (54),  using 

in , n 

(61)  and  by  applying  the  orthogonality  condition  (45): 


U - (n-1)  D 4'  . - (n+2)  D i|>  ..  , 

m,n  ra,n  m,n-l  in, n+1  in, n+1 


V - 1m  ill 
in,n  m,n 


(62) 


Note  that  the  expansion  for  1)  as  given  in  (60)  must  extend  one  degree 

above  that  defined  for  il'  . since  nonzero  values  of  Ui  i i i , are 

I in  | , | m 1 +M+ 1 

Implied  by  nonzero  values  of  'I'i  i i i • 

7 | m I , | m | +M 

Tlie  quantities  U , V and  V*tji  can  now  be  evaluated  at  points 
A ■ 2ttJ/N  , - arcain 
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where  J " 1 and  k “ • The  are  called  the  Gaussian 

latitudes.  Consider  for  example  V(X  ,ji  ) which  can  be  written: 

J * 


M . , r in  I +M 

ElmX  ' 

« n 


lm  P (u  ) 

m, n m,n  k 


m*  -M 


(6^) 


with  the  use  of  (41),  (46)  and  (62).  Similar  expressions  can  be 

written  for  U(A  and  V‘iji(Xj  .u^)  . The  outer  summation  can  be  carried 

out  very  efficiently  with  the  use  of  the  fast  Fourier  transform  method  which 

was  developed  by  Cooley  and  Tukey  (196S).  The  number  of  operations  required 

for  the  fast  Fourier  method  applied  over  N points  is  of  order  N Ion,  N , 

2 

while  for  the  direct  method  order  N operations  are  required.  The  fast 

Fourier  transform  method  Is  clearly  much  faster  than  tin  direct  method  for 

larger  values  of  N.  The  next  step  is  to  compute  UV2i|i  and  VV2\|i  at  each 

grid  point.  After  these  products  have  been  computed,  the  Fourier  transforms 

must  be  calculated  to  give  A and  R for  use  in  (66).  For  example, 

m m 

using  the  discrete  Fourier  transform: 


A (U  ) 
m k 


-imX 

e 1 <UV*t|»  . 


(64) 


where  -M  < m < M . A similar  expression  is  obtained  for  B (p,  ) . The 
fast  Fourier  transform  can  also  be  used  here  to  save  time. 

It  is  important  to  choose  N large  enough  to  avoid  aliasing  when  the 
products  are  transformed  back  to  wave  number  space  as  in  equation  (64). 
Orszag  (1969),  (1970)  suggested  that  N • 4M  would  be  needed,  but  later 
Orszag  (1971)  and  Machehauer  and  Rasmussen  (1972)  showed  that  N - 3M+1  was 
adequate  to  provide  alias-free  transforms. 
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can  be  computed  exactly 


Now  that  A (p.  ) and  B (p,  ) are  known,  F 

ra  k m k m,n 

from  (59)  If  the  degree  of  the  polynomials  is  less  than  or  equal  to  2K-1. 

The  maximum  degree  can  be  determined  from  (52)  by  noting  that  n is  a 

polynomial  of  degree  n and  by  considering  these  selection  rules  for  the 

interactions:  m * m^  + m^  , < 1^  + ‘ concll*sion  which 

is  given  in  Bourke  et  al.  (1977)  is  that  the  maximum  degree  is  5M-1,  so  that 

the  number  of  Gaussian  latitudes  is  K 5 M/2  . 

This  method  of  computing  F is  more  efficient  than  the  interaction 

m,n 

coefficient  method  and  it  requires  much  less  computer  storage.  The  number  of 

calculations  required  for  the  interaction  coefficient  method  is  of  order  (M^) 

3 

while  for  the  transform  method  it  is  of  order  (25  M ) (see  Bourke  et  al. 
(1977)].  It  will  be  shown  in  the  next  section  that  the  transform  method  is 
more  efficient  for  even  a moderate  value  of  M and  this  advantage  increases 
rapidly  with  M . 

7 Spectral  Model  of  Shallow  Water  Equations 

In  this  section  the  spectral  method  will  be  extended  to  the  primitive 
equations  and  it  will  be  demonstrated  that  semi-implicit  differencing  can  be 
applied  with  little  extra  effort.  The  shallow  water  equations  in  spherical 
coordinates  will  be  used  to  demonstrate  the  procedure  following  Eliasen  et  al. 
(1970r  and  Bourke  (1972).  The  equation  of  motion  and  the  continuity  equation 
can  be  written: 


-*  -*•  v*v 

« - (;+f)  k x v - v(<j>*  + Y)  * 

M’-  - v*$’v  - $6  . 


(65) 

(66) 


This  form  of  the  equation  of  motion  will  simplify  the  derivation  of  the  vor- 
ticity  and  divergence  equations.  Note  that  the  geopotential  has  been  split 

_ i 

into  a mean  $ , and  a departure  , which  will  facilitate  the  implementa- 

tion of  semi- implicit  time  differencing. 
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The  velocity  is  broken  into  rotational  and  divergent  parts  as  follows: 


V - k x V4>  + Vx  “ (U/cosqp)  i + (V/coscp)  T • 


(67) 


The  modified  components  U and  V will  also  be  used  here.  Now  form  the.  vor- 
ticity  and  divergence  equations  by  taking  7*  and  k*7x  of  (65)  which 
gives: 


3t 


- V-(c+f)  V , 


||  - k-Vxte+f)  v)  - v2  (4> ' + ^) 


(68) 

(69) 


The  vorticity  and  divergence  become 


C - V2*  , 6 =*  72X 


(70) 


In  spectral  models  it  is  convenient  to  replace  the  equation  of  motion  by 
the  vorticity  and  divergence  equations  because  the  relations  (70)  are  sim- 
plified when  spherical  harmonics  are  used  as  basis  functions.  This  form  of 
the  equations  is  also  more  convenient  for  application  of  semi-implicit  dif- 
ferencing. 

The  vorticity  equation  (68)  and  the  divergence  equation  (69)  can 
now  be  expanded  with  the  use  of  (67)  and  (70)  to  give: 


l ‘ - 


S"  (UV^)  + COSCP^  (VV2'I')1 


a cos  <p 


- 2fl(sin<p  V x + V/a)  , 


— V2v  . 

at  x 


V [£  (W2*)  - co.„£  (UV2*)] 


a cos  cp 


+ 2fi(sin<p7  - U/a)  - 7 ( 


2 2 
2 ,u  + V 


2 cos  <p 


+ ♦ ')  . 


(71) 


(72) 
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Similarly  Che  continuity  equation  (66)  becomes: 


If  w>  + cos*>  L W>1  - ^2X 


a cos  cp 


The  two  components  of  (67)  can  be  written: 

„ „ _ cogjg  M + I 

a 3<p  a 31  ’ 

v » I M + £°*SQ  IX 

a 3A  a 3qp 


(73) 

(74) 

(75) 


Equations  (71),  (72)  and  (73)  are  the  predictive  equations  for  \ 

and  and  (74)  and  (75)  are  diagnostic  expressions  for  U and  V . 

The  nonlinear  terms  in  these  equations  are  in  a convenient  form  for  the 
transform  method  which  was  presented  in  the  last  section,  since  the  multi- 
plication can  be  performed  at  the  grid  points  before  differentiation. 

Each  of  the  dependent  variables  is  expanded  in  terms  of  the  spherical 
harmonic  basis  functions  (41)  as  follows: 


( 76) 


( 78) 


These  expansions  are  for  the  rhomboidal  wave  number  truncation.  Equations 
(74)  and  (75)  are  transformed  in  the  same  manner  as  equations  (54)  were 
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In  the  last  section  and  the  result  is 


U - (n-l)D  ip  . -(n+2)D  ip  + imy 

m,n  m,n  m,n-l  m,n+l  m,n+l  m,n 


V — (n-l)D  Y . +(n+2)D  y + ii# 

m,n  ra,nAm,n-l  m,n+lAm,n+l  m,n 


Note  that  the  expansions  for  U and  V must  extend  one  degree  above  the 
expansions  for  ip  and  x * 

The  quantities  needed  for  the  nonlinear  terms  are  obtained  by  evaluat- 
ing the  sums  In  (76),  (77)  and  (78)  at  equally  spaced  points  in  longi- 

tude and  at  Gaussian  latitudes.  The  required  products  are  computed  at  each 
point  and  the  products  are  then  Fourier  transformed  in  longitude  as  follows: 


2\p  - a Am  elmA  , W2iji  - a Bm 


-“3£  c„'1"x-  v*’  -»3X)  D« 


2 2 
• »»■ 


U + V 2^  _ 

2 2^  E"> 


The  spectral  equations  are  formed  by  substituting  (76),  (77),  (78), 

(80),  (81)  and  (82)  into  the  system  (71)-  (73)  and  multiplying  each 

* 

equation  by  Y and  integrating  over  the  domain.  With  the  use  of  the 
ra,  n 

orthogonality  condition  (45)  the  equations  finally  reduce  to  the  following 


-n(n+l) 


3i|>m,n  1 I 1 

5t  'TJ  i-> 


r [imA  P - B 
2 m m,n  m dy 


+ 20[n(n-l)D  x , + (n+l)(n+2)D  .,X 

m,nAm,n-l  m,n+lAm,n+l 


- V_  _] 


-n(n+l) 
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m,n 

3t 


_ 1 
1-M 


, | ImB 


P 

m m,n 


dP 

+ A 

m dp 


dp-2fi[n(n-l)D  4/ 

m,n  m,n- 
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+ (n+1)  (n+2)D  .,4*  + 

ra,n+l  m,n+I 


V 1 

m,n 


f n(n+l)(E 


m,  ii 


> f ) 

m,n 
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n,n 


3t 


| imC  P 
m m,n 


dP 

D ) dU  + $n(n+l) 

m dll 


m,n 


(84) 

(85) 


where 


E 

m,n 


E 

-25 — P 

l-p2  m*n 


dp  . 


(«b) 


The  Integrals  are  evaluated  by  the  Gaussian  quadrature  formula  as  before, 
but  this  time  (5M+D/2  Gaussian  latitudes  are  required.  As  before  the 
required  number  of  longitudinal  grid  points  is  3M+1 . 

Bourkc  (1*172)  compared  the  efficiency  of  the  transform  method  to  the 
Interaction  coefficient  method  for  this  model.  Figure  2 shows  the  computer 
time  required  per  time  step  for  the  two  methods  as  a function  of  the  trunca- 
tion number  M.  The  figure  shows  clearly  that  even  for  M - 15  the  transform 
method  is  an  order  of  magnitude  faster  than  the  interaction  coefficient  method 
In  fact  the  Interaction  coefficient  method  becomes  almost  intractable  for  M 
much  larger  than  15.  At  M - 15  there  are  over  500,000  interaction  coeffi- 
cients. 

The  system  (83)-  (85)  is  very  convenient  for  the  application  of  serai- 

implicit  time  differencing.  All  terms  are  evaluated  explicitly  except  that 

$ in  (84)  and  y in  (85)  are  treated  implicitly.  These  two 
ih  « n ui  v n 

equations  are  easily  solved  for  (t+.\t),  and  equations  (83)  and  ( 84) 

ra,n 
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Figure 


2 

* -Computation  time  per  time  »tep  (»)  u a function  of 
cpectral  resolution.  Integrations  of  a global  spectral  modal 
employing  a transform  method  and  employing  the  interaction 
coefficient  method  are  oompared. 


can  then  be  solved  explicitly.  In  contrast  finite  difference  models  require 
the  solution  of  a Helmholtz  equation  for  $(t+At)  , at  every  time  step 
Thus  in  spectral  primitive  equation  models  a much  longer  time  step  can  be 
used  with  almost  the  same  computat ional  effort  per  time  step. 

The  introduction  of  the  transform  method  and  semi- implicit  differencing 
have  made  the  spectral  primitive  equation  models  competitive  with  finite 
difference  models  for  global  prediction.  _ The  procedures  used  in  this  section 
are  easily  extended  to  baroclinic  models  as  has  been  done  by  Bourke  et  al. 
(1977),  Machenhauer  and  Daley  (1972)  and  Hoskins  and  Simmons  (1975).  Com- 
parisons have  shown  that  as  good  or  better  forecasts  can  be  made  with  global 
spectral  models  than  with  finite  difference  models  which  use  the  same  amount 
of  computer  time  (Doron  et  al.  (1974)  and  Daley,  Girard  and  Simmonds  (1976). 

It  should  be  pointed  out  that  energy  is  not  exactly  conserved  in  this 

model  even  with  continuous  time  variation.  This  is  because  the  kinetic  energy 

*♦ 

for  the  shallow  water  equations  is  proportional  to  <J>V • V which  is  a cubic 
energy  form,  and  consequently  the  analysis  of  Section  3 does  not  apply. 
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However  the  nonlinear  terms  are  computed  very  accurately  in  spectral  models 
and  experience  shows  that  the  energy  is  in  fact  very  nearly  conserved. 
Bourke  (1972)  integrated  the  model  which  was  developed  in  the  section  for 
116  days,  and  obtained  an  energy  change  of  only  2 percent. 

8 Advection  Equation  with  Finite  Elements 

In  this  section  the  finite  element  method  with  linear  elements  will  be 
applied  to  the  advection  equation 


3u 

3t 


+ c 


3u 

3x 


(87) 


This  equation  has  been  treated  extensively  with  various  finite  differ-' 

v 3 

ence  schemes.  The  Calerkin  equation  is  obtained  by  setting  ■ c ^ in 
( 22)  which  gives 


*1 


dx 


0 


1-1.. 


,N. 


(88) 


The  linear  basis  functions  cpj  (x)  are  defined  by  (12)  and  a typical  one 
is  shown  in  Fig.  1.  In  this  application  u is  periodic  so  that  the  basis 
functions  must  satisfy  cpo  - and  cp.  * ^+1  • 

The  first  term  in  (88)  can  be  evaluated  from  (17)  which  is  of  the 
same  form,  and  the  second  term  can  be  computed  with  the  use  of  (14).  The 
resulting  equation  with  i ■ m can  be  written: 


, du  . , du  du 

1 / — + * _J»  + _J 

6 vdt  dt  dt 


m-!)  Vl  " Vl  . o 

' 2Ax 


(89) 


The  advection  term  is  the  same  as  is  obtained  from  centered  differencing, 
but  the  time  derivative  appears  as  a weighted  average  over  three  points.  It 


will  be  seen  later  that  this  greatly  increases  the  accuracy  of  the  solution. 


Now  apply  leapfrog  time  differencing  which  gives  the  following  equation: 


;(u 


.+4(u  _.,-u .)+u 


)+  -^(u  _-u_  , ) - 0 


12Atv  ra+l,n+l  nrfi,n-l  m,n+l  ra,n-l  m-l,n+l  2Ax  nrfltn  m-lfn 


(90) 


The  stability  and  phase  error  can  be  investigated  by  substituting  u 
A exp[i(yAxra  + oAtn) ] into  (90).  This  yields 


m,n 


sinoAt  ■ -(cAt/Ax)(3  sinyAx)/(2  + cosyAx) 


(91) 


The  solution  is  stable  (i.e.  (sinoAt)  1)  if 


3cAt/Ax[sinyAx/(2  + cosyAx)]  < 1 . 

max  — 


The  bracketed  term  is  a maximum  when  yAx  =*  120®  , so  that  the  stability  con- 
dition becomes 

|cAt/Ax|  < 1//3  . (92) 

This  criterion  is  considerably  more  restrictive  than  the  CFL  condition  which 

arises  from  the  leapfrog  finite  difference  scheme.  However  it  will  be  shown 

that  (90)  gives  even  better  phase  speed  than  the  fourth-order  leapfrog 

scheme  for  which  the  computational  stability  criterion  is  ]cAt/Ax|  <^0.73 

Thus  it  is  not  unreasonable  that  -the  leapfrog  finite  element  scheme  would 

have  a more  restrictive  computational  stability  criteria. 

The  finite  element  formula  with  leapfrog  time  differencing  is  actually 

implicit,  since  the  new  value  u ,,  cannot  be  obtained  explicitly  from  the 
r ra,n+l 

earlier  time  values.  Thus  it  seems  reasonable  to  use  a fully  implicit  form 
which  does  not  have  the  timestep  restriction  (92).  Consider  the  following 
time  difference  approximation  to  (89): 
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(93) 


:(u 


.-u 


+ 4(u ,,-u ) + u 


■»  .1  -> 


6At  nrt-l,n+l  nrH,n  m,n+l  ra,n  m-l,n+l  m-l,n 


c 

+ 77—  (u  ,,-U  , . + u -u  , ) ■ 0 

4Ax  m+l,n+l  m-l,n+l  m+-l,n  m-l,n 


This  fully  Implicit  scheme  can  be  shown  to  be  neutral  for  all  timesteps, 
and  it  should  require  about  the  same  effort  per  time  step  as  (90).  For 
this  reason  implicit  time  differencing  schemes  are  often  desirable  when 
finite  elements  are  used. 

The  phase  speed  for  the  leapfrog  scheme  is  given  by 


, 1 . cAt  3 sinyAx. 

^ " - a/u  = Wt  arcsln  ~ Ax  T+ToiW  • 


(94) 


If  At/Ax  and  p are  fixed,  this  expression  approaches  c as  At  -*•  0 , 
which  shows  that  the  solution  converges.  If  At  is  small  in  comparison  to 
Ax/c  , this  formula  reduces  to 


c . c 3 sinyAx  m c 
F uAx  2+cospAx  yAx 


sinyAx 


[1-2/3  sin  (yAx/2)] 


(95) 


Table  1 contains  c/c  from  ( (95)  for  typical  values  of  L 

F 


L 

2Ax 

3Ax 

4Ax 

6Ax 

FEM 

0 

0.83 

0.96 

0.99 

4th  order 

0 

0.61 

0.85 

0.96 

Table  Is 

c/cv  for 

the  FEM 

solution 

and  for  4th  order 

differenced  scheme  for  various  wavelengths  L. 


The  table  also  includes  the  ratio  for  the  fourth  order  scheme  from  the  limit 
for  small  At  . The  finite  element  formula  (95)  can  be  expanded 
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in  pAx  which  leads  to  an  error  that  is  of  order  (pAx)  . Table  1 shows 
that  although  the  linear  finite  element  equation  and  the  fourth  order  finite 
difference  equation  have  the  same  order  of  truncation  error,  the  finite 
element  equation  is  much  more  accurate.  At  L “ 3Ax  the  finite  element 
solution  gives  only  17Z  error  in  phase  speed,  while  the  fourth-order  finite 
difference  gives  393.  However  for  L • 2Ax  , c - 0 , which  indicates  that 
the  finite  element  computational  group  velocity  is  very  large  for  this  wave- 
length. This  can  be  shown  by  differentiating  as  follows: 


c 


d (ucp) 

dp 


[ 2cospAx  + 1 ] 

2 

(cospAx  + 2) 


( 96) 


When  L - 2Ax  (pAx  - tt)  this  formula  gives  G - -3c  which  is  much  larger 
than  the  -(5/3) c that  occurs  with  fourth-order  differencing.  This  suggests 
that  small  scale  noise  will  propagate  very  rapidly  in  finite  element  models. 

This  tendency  toward  noisiness  has  been  observed  in  various  finite  element 
models.  The  degree  of  accuracy  indicated  above  for  the  finite  element  model  has 
been  verified  by  Cullen  (1973)  in  a two-dimensional  advective  problem.  It 
should  be  noted  that  although  the  FEM  gives  a solution  for  all  values  of  x 
in  the  range  considered,  the  high  accuracy  is  only  obtained  at  the  nodal 
points  since  the  fields  are  assumed  to  be  linear  between  nodal  points.  In 
the  next  section  the  method  will  be  applied  to  the  barotropic  vorticity  equa- 
tion. 


39 


9 . Barotropic  Vorticity  Equation  with  Finite  Elements 


In  this  section  the  finite  element  method  will  be  applied  to  the  non- 
linear barotropic  vorticity  eqviation  in  a two-dimensional  domain.  The  basis 
functions  will  be  linear  functions  on  triangular  elements.  The  barotropic 
vorticity  equation  can  be  written 


In 

3t 


-k  x Vt^*Vri  , 


97) 


where 


n - f(y)  + V2^  , 


(v  98) 


is  the  absolute  vorticity. 

Following  Fix  (1975)  both  tji  and  n are  expanded  in  terms  of  the  basis 
functions  cpj(x.y)  as  given  below: 


N 


iKx 


,y,t)  = >IVj(t)  cpj(x.y)  , 


j-l 

N 


n(x,y,t)  =■  Hj(t)  cpj(x,y) 

J = 1 


( 99) 


( 100) 


When  the  Galerkin  method  is  applied  to  ( 98)  the  following  is  obtained: 

^ (t)  f f CpiV2CPj  dA  = - f f tp1(x,y)f  (y)dA  + ^ TTj  (t)  I I cpi<x,y)  cpj  (x,y)dA  , 
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< 

i 
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for  1*1,.., N.  Since  linear  basis  functions  will  be  used  it  is  necessary  to 
integrate  the  left  hand  side  by  parts,  which  gives: 


The  boundary  terms  which  arise  from  the  integration  by  parts  were  set  to  zero, 

by  assuming  that  either  ip  is  periodic  in  space  or  that  there  is  no  flow 
• “►  — ► ■+■ 

normal  to  the  boundaries(i.e.  ,kxVi(;*n  = 0,  where  n is  a unit  vector  normal  to 
the  boundary).  Now  apply  the  Galerkin  method  to  the  vorticity  equation  ( 97), 

which  leads  to  the  following  form: 


Z)  sr1/  Ky*  ■ - X)  2 Vk jjv !**%,•  'va  • 


for  i=l,...,N.  This  equation  is  of  the  same  form  that  was  obtained  with  the  spectral 
model,  but  the  nonlinear  term  requires  much  less  effort  because  the  only  cp's 
which  interact  are  those  which  are  physically  adjacent. 

The  equations  ( 101)  and  ( 102)  conserve  both  mean  square  vorticity 

(enstrophy)  and  kinetic  energy.  The  enstrophy  conservation  can  be  shown  by 
multiplying  ( 102)  by  and  summing  over  i . When  the  sumations  are 

taken  under  the  integrals,  the  form  ( 26)  is  founo  Since  the  integral  of 

Hk  x Vtp^Vri  vanishes,  the  conservation  of  0^/2  follows  directly.  The 
kinetic  energy  change  can  be  examined  by  first  differentiating  ( 101)  and 

substituting  the  result  into  ( 102)  which  gives: 


Vcp^’VcpjdA  = - 


Vcf^dA  , 


( 103) 


for  i=l,2,...,N.  Multiply  this  equation  by  and  sum  over  i . The 

resulting  equation  is  again  of  the  same  form  as  ( 26)  and  the  left  hand  side 
is  the  derivative  of  the  total  kinetic  energy.  Since  the  integral  of 
tpk  x ViJi*Vn  is  zero,  the  energy  is  conserved.  These  results  are  not  dependent 
on  the  particular  basis  functions  which  are  employed. 

The  systems  of  equations  ( 101)  and  ( 102)  can  be  written  in  matrix 

forms  which  are  more  convenient  for  solution.  Let  iji  and  ri  be  column 


k i 
' . 


. 


i 


1 


\ 

i 


^ i 


vectors  of  the  values  of  4<^  and  rt^  , respectively.  Then  ( 101)  takes 

the  form: 

K< p m Q , (104) 

where  the  elements  of  the  matrix  K arc 

Kt)  “ ff  Vf4J  * VqVU  ' ( 105) 

•> 

and  Q Is  a column  vector  of  the  right  hand  side  of  ( 101).  Similarly, 

system  ( 102)  becomes 

M ^ - J , ( 106) 

where  the  elements  of  M are 

M - ff  cp(Cpj  dA  , ( 107) 

and  .1  is  a column  vector  of  the  right  hand  side  of  ( 102). 

The  solution  procedure  will  be  illustrated  for  the  case  where  leapfrog 
time  differencing  Is  used  In  ( 105)  which  leads  to  the  equation: 

MAn  - 2At  J , (108) 

> ► > 

where  An  = n , - n , . The  matrices  K and  M are  computed  Initially 
n f 1 n- 1 

and  stored  for  later  use.  The  equations  can  be  Integrated  beginning  with 

4'  , ti  and  rj.  • The  right  hand  side  of  ( 108)  can  be  computed 

J , n J , n- 1 J.n 

from  4i  and  n , and  that  equation  can  then  be  solved  for  An.  . This 
J.n  J.n  J 

Increment  can  be  added  to  p,  , to  obtain  rt  . , • With  these  values  the 

J,n-1  J .n+1 

right  hand  side  of  ( 104)  can  be  computed,  and  ( 104)  can  be  solved  for 

4'  , and  the  process  can  be  continued.  In  this  procedure  It  Is  necessary 

J.n  + l 

to  invert  the  matrices  K and  M during  each  time  step.  These  matrices  are 
very  sparse  since  only  adjacent  elements  interact.  In  some  cases  direct 
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methods  can  be  used,  but  iterative  methods  are  much  more  flexible^ 

Cullen  (1973)  has  shown  that  the  two  dimensional  advective  stability 
criterion  for  linear  elements  Is 


( 109) 


where  d Is  the  distance  between  nodal  points.  This  Is  consistent  with  the 
one-dimensional  result  (92),  because  the  step  from  one  to  two  dimensions 
Is  usually  achieved  by  replacing  the  grid  si-^e  with  d//2  . Tn  this  applica- 
tion [ c|  would  correspond  to  the  maximum  velocity  In  the  domain.  Since  the 
condition  (108)  is  rather  restrictive  for  At  and  since  two  matrices  must 
be  inverted  per  time  step  it  may  be  worth  while  to  use  a fully  Implicit  form 
similar  to  (93). 

The  natural  generalization  of  the  tent  function  in  one  dimension  to  two 
dimensions  is  a basis  function  which  is  composed  of  triangular  elements.  On 
each  triangle  the  function  varies  linearly  from  0 at  two  vertices  to  l at  the 
third  which  is  the  nodal  point  for  the  basis  function.  Figure  3 shows  how 
a typical  basis  function  cpj  is  constructed  on  a rectangular  grid  of  nodal 
points.  This  function  is  the  sum  of  the  six  plane  surfaces  that  are  associated 
with  each  triangle.  The  basis  functions  can  be  equally  well  constructed  when 
the  nodal  points  are  irregularly  located,  and  it  is  not  necessary  to  have  six 
triangular  elements  In  the  construction. 

The  elements  in  the  matrix  equations  (104)  and  (106)  are  obtained  by 
evaluating  the  integrals  in  equations  (101)  and  (106).  These  integrals 
can  be  reduced  to  a series  of  Integrals  over  triangles  such  as  are  shown  in 
Figure  3.  Within  each  triangle  any  point  is  affected  by  only  the  three 
basis  functions  which  have  nodal  points  at  the  three  vertices  of  the  triangle. 
Zlcnkiewlcz  (1971)  and  Desal  and  Abel  (1972)  describe  a convenient  procedure 


for  evaluating  the  integrals  over  each  triangle.  This  Involves  introducing 
triangular  coordinates  which  vary  linearly  across  each  triangle  in  the  same 
manner  as  the  basis  functions.  The  integrals  can  then  be  evaluated  quite 
generally. 

A.  rigorous  mathematical  analysis  of  the  finite  element  method  is  given 
in  the  book  by  Strang  and  Fix  (1973).  The  stability  and  convergence  of  the 
method  are  discussed  in  considerable  detail.  Most  finite  element  applications 
are  based  on  a variational  formulation  rather  than  the  Galerkin  approach  which 
has  been  used  here,  although  the  Galerkin  method  is  most  appropriate  when  time 
dependence  is  included.  Finder  and  Gray  (1977)  developed  the  finite  element 
metood  with  the  Galerkin  approach,  and  gave  applications  in  hydrology  which 
has  similar  equations  to  those  which  occur  in  numerical  weather  prediction. 

The  finite  element  method  has  been  applied  to  atmospheric  prediction  with 
the  primitive  equations  in  shallow  water  form.  Cullen  (1974)  and  Hlnsman  (1975) 


Fig. 
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Construction  of  the  basis  function 
rectangular  array  of  nodal  points. 


on  a 


carried  out  global  forecasts  with  these  equations  using  linear  basis  functions 
on  triangles  as  discussed  in  this  section.  The  elements  were  efficiently 
arranged  so  that  the  area  of  each  element  was  almost  the  same  over  different 
parts  of  the  globe.  Most  global  finite  difference  models  have  a large  varia- 
tion in  grid  size  between  the  equator  and  the  pole,  and  consequently  are  not 
very  efficient. 
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Stanlforth  and  Mitchell  (1977)  reformulated  the  shallow  water  equations 
In  terms  of  the  vorticity  and  divergence  as  was  done  In  section  6.6  for  the 
spectral  model.  In  this  form  semi-implicit  time  differencing  can  be  applied 
easily,  which  allows  a much  larger  time  step.  This  is  very  important  since 
the  finite  element  method  generally  requires  more  computer  time  per  time  step. 
Stanlforth  and  Mitchell  also  found  very  little  noise  generation  in  their  fore- 
casts, whereas  many  finite  element  primitive  equation  models  tend  to  generate 
small  scale  noise  if  no  smoothing  is  used  [Cullen  (197i>]. 

The  finite  element  method  when  applied  to  meteorological  equations  gives 
very  accurate  phase  propagation  and  also  handles  nonlinearities  very  well. 

The  main  drawback  to  the  use  of  the  method  is  the  requirement  that  an  equation 
solver  must  be  applied  to  invert  a large  matrix  at  every  time  step  for  every 
variable.  The  development  of  flexible  exact  solvers  for  these  matrices  is  of 
great  importance.  The  finite  element  method  can  easily  be  applied  to  variable 
resolution  problems,  but  some  finite  element  models  do  tend  to  produce  noise 
probably  as  a result  of  the  large  spurious  group  velocity  for  the  shortest 
wave.  However,  the  formulation  of  Stanlforth  and  Mitchell  (1977)  seems  to 
reduce  this  problem  considerably.  Schoenstadt  (1978)  has  shown  that  noise  is 
generated  in  finite  element  models  where  all  variables  are  carried  at  the  same 
models  points.  When  the  variables  are  staggered  at  different  model  points  or 
when  the  vorticity  and  divergence  equation  are  used  this  problem  can  be 
avoided.  The  general  procedure  used  by  Stanlforth  and  Mitchell  (1977)  appears 
to  be  superior  because  semi- implicit  differencing  can  be  easily  implimented, 
and  the  forecasts  are  not  noisy. 
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